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Abstract 

An efficient new method is presented to calculate the quantum transports using periodic bound- 
ary conditions. This new method is based on a method we developed previously, but with an 
essential change in solving the Schrodinger's equation. As a result of this change, the scattering 
states can be solved at any given energy. Compared to the previous method, the current method is 
faster and numerically more stable. The total computational time of the current method is similar 
to a conventional ground state calculation. Details of the procedure is presented in the current 
paper. 

PACS numbers: 71.15.-m, 73.63.-b, 73.22.-f 



I. INTRODUCTION 



Molecular elastic quantum transport has been studied intensely in recent years both in 



theory and experiment 



mmm 



Theoretically, the total current through a molecule 



connected by two electrodes can be calculated as: 

2e ff^R 
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Y.Tn{E)dE, (1) 



where /i^ and fiR are left and right electrode Fermi energies (assuming the current flows 
from right to left in z direction), and Tn{E) is the transmission coefficient for the nth right 
hand electrode band at energy E. One common way to calculate Tn{E) is to calculate the 
scattering states ipsd^) at energy E which satisfies the Schrodinger's equation: 

Hij^r) = E^sc{r) (2) 
while having the following boundary conditions: 
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Note that, in Eq(2), H = { — + V{r) + Vnonioc} is the single particle Hamiltonian. In 
Eq(3), (f)^^^\r) = Un,k„{''^)exp{ik^^^^ z) , are the right going running waves in the the right(R) 
and left(L) electrodes, and (f)n^^^* are the left going running waves. E^^^^k^^^^) = E are 
the electrode band structure. The summation En in Eq(3) stands for all band n and k^'^^^ 
which satisfy E^'^^\k^'^^'>) = E. In Eq(3), we have assumed that dE^^^\k)/dk > for band 
n. If one band n has dE^^^\k)/dk < 0, then its corresponding (p^^^^r) in Eq(3) should 
be replaced by 0^^^^*(r) for that particular band. Eq(3) describes an incoming running 
wave 0^*(r) from the right electrode band m which is scattered back through outgoing 
running waves B^(f)fl^{r) at the right electrode, and transmitted into the left going running 
waves A^0^*(r) at the left electrode. As a result of this scattering process, the transmission 
coefficient for channel m can be calculated as 

TUE) = E \A^\'{dE^{k)/dk)\,=,r,]/{dE^{k)/dk)\,=,n. (4) 

n 

Note that, due to current conservation, we have: 



Y: \A^\'{dE^{k)/dk)\,=,. + Y: \B^\\dE^{k)/dk)\,=,n = {dE^{k)/dk)\,=,n. (5) 

Normally, Eqs(2),(3) are solved by using transfer matrix method or Lippmann- 

Schwinger equation (5|. However, transfer matrix method could be unstable in a multi- 
channel electrode and it is difficult to deal with the nonlocal pseudopotentialjsl, and the 
application of Lippmann-Schwinger equation is computationally expensive^]. 

In a previous publication j^, we have described a new method to calculate Eqs(2),(3) 
which uses a supercell with periodic boundary conditions, just like in a conventional ground 
state total energy calculation. In that method, a supercell eigenstates are solved using con- 
ventional conjugate gradient methods Q]. Then perturbations at one end of the electrode are 
introduced, and the eigenstates are recalculated using the same conjugate gradient method 
j^. Next, these eigenstates are linearly recombined to make it satisfy the boundary condition 
Eq.(3). Although highly efficient compared to methods before it, and simple to implement 
since it uses only conventional ground state codes, that method has a drawback. In that 
method, the energy E in Eq(2) can only be the eigen energies of the original supercell 
Hamiltonian H. As a result, Tn{E) is only known for a finite number of E (or say To 
overcome this problem, one needs to fit Tn{kn) with a continuous function before it is used 
to calculate the total current in Eq(l). Although it has been shown in Ref. that this 
fitting over the fc„ points is equivalent to the k-point summation in a supercell ground state 
calculation and it works fine in the case considered, but there might be cases where denser 
energy E points are needed, for example, close to a weakly coupled resonant tunneling. In 
this paper, we will provide an essential modification over our previous method. Under this 
new method, Eqs(2),(3) can be solved for arbitrary E, and the overall computation is faster 
than the previous method for much denser E point grid. We also provide details of the 
whole procedure in the current paper. 

II. THE FORMALISM 

In order to compare with our previous method in Ref. 0, we choose the same svstem as 
studied in that paper. The system is schematically shown in Fig.l taken from Ref. In the 
system, a benzene molecule is connected by two Cu quantum wires through the bonds of two 



sulfur atoms. In Fig.l, left hand side B and right hand side B are periodically connected. 
As in Ref. 0, we will do nonselfconsistent calculations for finite bias, although selfconsistent 
calculation is straight forward using the scattering state solutions of Eqs(2),(3). To get the 
potential of a finite V bias system, self-consistent local density approximation (LDA) ground 
state calculation is performed for the zero bias system. Then, an additional smooth function 
is added to the potential to raise the right hand electrode by V/2 and lower the left hand 
electrode by —Vjl. Norm conserving pseudopotentials are used, so is a planewave basis set 
with a 30 Ryd cutoff. 

With some minor changes of notations, the essential idea of Ref. [its Eq(6)] can be 
recast as to solve the wavefunction '^[[){r^ : 

(if-E)^(0(r) = W^(0(r), (6) 

here are some perturbation functions which are only nonzero away from the 

molecule as shown in Fig.l. Note that, 4'{i){r) has a supercell k-point K^, thus = 
ip(i){r)exp{—i'Kzz) is periodic. could be, for example 7i/2Lz where is the supercell 
length in the z direction. After a few ip{i){''^) are solved for a same energy E, these ip{i){r) can 
be linearly recombined with proper coefficients to generate a scattering state ipsc{f^) which 
satisfies the boundary conditions of Eq(3). In Ref. |8|, Eq(6} is obtained by combining two 
eigen state equations with two different Wm [ Eq(6) in Ref. 13] • Its advantage is that it needs 
only conventional eigen state calculations, thus there is no need to change a ground state 
code. The disadvantage, however, is that the energy E in Eq(6) can only be the eigen state 
energy Ei of the unperturbed = 0) supercell system. Here, we will solve the Eq(6) 

directly using the conjugate gradient method. The approach is very similar to the method 
used in perturbation linear response theory 

Notice that, the linear equation (6) can be rewritten as an optimization of the following 

F: 

F =< ^/j^i)\H - > - < ^/Jii)\W(i) >-< W^dI^^d > . (7) 

Preconditioned conjugate gradient method can be used to solve the minimum of F. Un- 
fortunately, for an arbitrary E, the matrix H — E in the above equation is not positive 
definite, which makes the conjugate gradient method diverges. However, this problem can 



be circumvented if we have the eigenstates {ipi, Ei} of the unperturbed Hamiltonian H. 
First, if we know all the eigenstates {ipi,Ei} of H, the in Eq(6) can be solved directly 
as: 



i ^ 

In practice, however, we usually only solve the eigenstates {ipi, E^} (using the conjugate 
gradient method P]) up to an energy E' (with E' > hr). Let's denote this converged set of 
eigenstates as i = 1, A^. Then the idea is to deflate these N eigenstates from both ip(j_-^ and 
W{i), and solve the remaining part of the wavefunction. More specifically, we can define: 

N 

^(0 = ^^(0 = ^(0 - < V^i|^(0 > V'i (9) 

i=l 

and 

AT 

W[i^ = = Wil) - ^ < iP,\W^i) > ij,. (10) 

1=1 

Then the linear equation 

(if-E)^f)(r) = iy(f)(r), (11) 
in the subspace of projector P can be solved as the minimum of 



=< ^f,^\H - > - < ^P^\wi;^ > - < ty^f^l^f) > . (12) 

Note that, now the effective matrix P{H — E)P is positive definite as long as E is 
lower than E^. When using conjugate gradient method to solve Eq(12), the projector 
P are repeatedly applied to the wavefunctions and search directions, so that the whole 
minimization is done within the subspace defined by P. After ip^^-^ is solved, of Eq(6) 
can be obtained from Eqs(9) and (8) as: 

«...-«f, + E^S^*. (13) 

The convergence of ip^^-^ under the conjugate gradient method is very fast since the effective 
band gap for ip^^-^ is E]^ — E. Kinetic energy G-space diagonal preconditioning can be used, 
just as in the conventional ground state conjugate gradient method^. Fig.2 shows a typical 
convergence for a ip^i^^ state. It shows that 20 conjugate gradient line minimizations is enough 



to converge the V'^) to 10~^ (a.u.) starting from zero. The ip^i) can be converged to the same 
accuracy as that of {ipi, Ei}. We find similar convergence for all and E. 

To linearly combine to generate ipse of Eq(3), we want to have linearly independent 
?/'(;). According to Eq(8), if M W(^i-j are linearly independent, then the corresponding M ipii^ 
are also linearly independent. As in Ref. jSj, we intend to use the V point electrode states 
as W{i). However, if one or few Ei are very close to E, then their corresponding ipi terms 
might dominate the expression in Eq(8), which can make the ipij^) lie in similar directions 
for different /. To avoid this situation, we have modified W{i) as following. Let's use to 
denote the F point electrode states in real space and here / is an index for different bands. 

are only nonzero at the last primary cell of the right electrode as shown in Fig.l. Then 
from we can generate W(i) using the following iterations from = 1 to = / — 1: 

- H/(/),^ - l^M— 7 TTTr — (14j 

< Wm^ I > 

and W(^i)^i = W^i^ and W(^i) = W(^i)^i. In the above equation, is the i which gives the 
maximum | < ipi\W(^^) > /{Ei — E) \ for a given /i. Note that, it is easy to show from Eq(14) 
that < 4'mjW(^i) >= for all the < 1. In other words, ■?/'(/) as described in Eq(8) [or 
Eq(13)] will not have the ipi component if tpi is a maximum component in one ^/'(^) with 
< /. This makes and ■?/'(^) (/i < /) unlikely to lie in very close directions. Note that 
W(^i-j from Eq(14) will still only be nonzero in the last primary cell of the right electrode. 

After the M ipii) of Eq(6) are calculated following the above procedure for a given E, 
we will combine these to generate the scattering states ipse of Eqs(2),(3). This part is 
similar to what we have done in Ref. jsi However, more details will be provided here. A 
band structure alignment between the right and left electrodes is illustrated in Fig. 3 with an 
4 V bias. The numbers in the right electrode band structure are the index of the electrode 
bands. In our calculation, we have over cautiously used 11 bands as W(i) in Eq(6). As a 
result, for a given energy E, we will have 11 ipi^iy Since ip*^i-^ (which has a — instead of K^) 
also satisfy Eq(6), we end up having 22 wavefunctions to be used in the linear combination 
to generate ipse [in the following, we will denote all these 22 wavefunctions as ip{i), with 
/ = 1,M and M being 22]. As described in Ref. [sl, we will first decompose each ip(i) at 
left and right electrode primary cells Vt^ and fi/j by the electrode wavefunctions. As shown 
in Fig. 3, for a given energy E, we can find the corresponding k^{E) and k^{E) (the small 



black dots in Fig.3 on the dashed hne E). If the numbers of k in the left and right electrodes 
are Nl and Nn respectively, then there will be 2{Nl + Nji) electrode running waves states 
(counting also the —k). Like in Eq(3), let's use (f)^^^^ and ipn^^^* denote these electrode 
running wave states, then the decomposition of can be written as: 

, Ji::=fiK(OeW + 5„^(O0^(r)] iizenn 

ipmir) = < (15) 

[ E:ii[A^(O0^*(r) + B^^iimr)] ifzen^ 
The coefficients A^, can be calculated by the overlap matrix < > and 

the projection matrix < ip(^i)\(l)^^*^ >. The same for Af^, B^. The electrode wavefunctions 
0„(r) are pre-calculated at 50 k points. Then the 0^(r) for a given k^ point is obtained vs 
interpolation. 

Now, combining 'ip(i){r) we have the scattering state as: 

M 

^.c(r) = E^^V^(0W- (16) 
1=1 

Note, ipsc{f) satisfies the Schrodinger equation (2) in the region where W^(/)(r) are zero 
(away from B in Fig.l). According to Eq(15), we have ipsci'f') at and Ql as: 

ipscir) = \ (17) 
E^;Si{Efii A^(/)Q]0^*(r) + lEfiiBf^mmr)} if ^ e fi. 

Comparing this equation with the boundary equation (3), we have the following Nj^ + N^ 
linear equations for a scattering state based on a 0^* incoming wave: 

Efii A^(OG = 5„,^ forn = l,iV^ ^^^^ 

EfiiB^:il)Ci = forn = l,iV^ 
Note that if M > Nji + there can be a solution for the above equation. Given the 11 
bands we used as W(^i) (M = 22), we find that this is always true. When M > Nji + Nl, 
Eq(18) is under determined, meaning there are more than one solutions of Ci. In this 
case, it makes sense to require the minimum of Efii \Ci? while Eq(18) is satisfied. This 
linear algebra problem can be solved using standard numerical routines, like the ZGELSS 
in LAPACk|io|. 

After Eq(18) is solved, then we will have a scattering wave ipse, which satisfies the 
Schrodinger's equation (2) within the region from to Qr, and the boundary condition 



of Eq(3) at Ql and Qr. We can discard the tpsc of Eq(16) for the regions outside Ql and 
Qr (near boundary B). Instead, the real scattering state can be extended following the 
propagations of the nonzero electrode running waves in Eq(17) into negative and positive 
infinities. As a result, the boundary conditions of Eqs(17),(18) at Qr and Ql are the same 
as the boundary conditions of Eq(3) at 2; — > oo and z — > —00. 

III. THE EVANESCENT STATES 

Above discussions are complete if the electrodes are sufficiently long, so there are no 
evanescent states at and D,r. In practices, however, we found evanescent states often 
exist. There are two type of evanescent states. The first type (type I) originates from the 
molecule and decays out in the electrodes (i.e. e~'^^ in the right electrode and e**^ in the left 
electrode). The second type (type II) originates from the artificial boundary B in Fig.l, and 
decays towards the molecule (i.e. e**^ in the right electrode and e~'^^ in the left electrode). 
While the first type evanescent states could be physical, existing in a scattering state ipse, 
the second type of evanescent states are artificial due to our use of boundary condition and 
perturbation W(^i) near B . If we can calculate the running wave coefficients A^'^^^l) and 
B^^^\l), then even if we ignore the evanescent states in Eq(15), our resulting scattering state 
Ipse constructed from Eq(16) and Eq(18) will still be correct. This is because the evanescent 
states can be added in as additional terms in Eq(17). As we extend our boundary condition 
from and D,r to —00 and 00, the first type evanescent states will decay out, and we 
can simply remove (subtract out) the second type evanescent states without affecting the 
Schrodinger's Equation (2) (assuming its amplitude near the molecule is sufficiently small). 
As a result, we will still have a boundary condition as in Eq(3). 

However, it is helpful to include the evanescent states in the decomposition Eq(15) for 
two reasons: (1) To accurately calculate the running wave coefficients A^^^\l) and B^'^^\l)] 
(2) In Eq(16), to avoid the case where large artificial second type evanescent states exist 
and dominate the equation. As a result, they have significant tail amplitudes near the 
molecule (compared to the running wave amplitudes) . If this is true, then these second type 
evanescent states cannot be simply removed without introducing errors. 

To get A^{1) and B^{1) from Eq(15), we have used the overlap matrix < (t>n\4'm ^^ri 
< (t>n*\(t>m >nR: < (l>n*\(t>m >^R ^nd projections < iJi\(t)^ >nj, and < >nR: then 



solved the resulting linear equation. Here the subscript means the integration is done 
only within Q^. Since the running waves and evanescent states are not orthogonal within 
Qr, then ignoring the evanescent states in Eq(15) will introduce errors in the resulting A^[l) 
and B^{1). The situation is the same for the left electrode. 

The evanescent states are originated from the real kg points when dEn{k) / dk\k=k^ = 0. 
Here the subscript e stands for evanescent states. This happens at the F and X' points of 
the electrode band structure as shown in Fig.3, and at the place where two bands anticross 
each other and form a band gap. At the F point, an evanescent state line runs downward 
starting from a real F point band structure energy. At the X' point and other anticrossing 
points, two evanescent state lines connect the two real k^, points at the opposite edges of 
an e nerg y gap. A detail description of the complex band structure is given by Chang in 
Ref. I3- Although calculating the complex band structure is possible it is difficult 
for a nonlocal pseudopotential Hamiltonian as is used here. As a result, we have used the 
real k^, point Bloch states (f)n,e{r) = 'iJ-n,ke^xp{ikez) to approximate the evanescent states. 
Notice that, these are just normal running wave states, except that they carry no current 
since dEn{k) / dk\k=k^ = 0. A more accurate approximation is to add an exponential decaying 
factor exp{K,z) or exp{—Kz) to 0„ e('")- However, since we are only going to use (f)n,e{r) within 
Qn OT Ql, and the surviving evanescent states within or Ql should have a small k, we 
found it is okay for not using these decaying factors. By not adding this decaying factors, 
we also do not distinguish the type one and type two evanescent states. 

Unlike the running wave states whose number is finite for a given energy E, there can be 
many (actually infinite if we have an infinite basis set) evanescent states for a given E. This 
is because at the F point of the band structure, every new band will have an evanescent 



state line running downward in energy 



12j |. As a result, in Eq(15), we cannot include all 



the possible evanescent states 4>n,e{r) = Un^k^expiikez) for a given E. On the other hand, 
in practice, it is not necessary to include the evanescent states which are originated from 
running wave energies En{ke) which are far away from E, because they will have fast decay 
factors exp{K,z), thus should not exist in Ql or fi/j. Because of this, we have the following 
practical procedure in solving Eq(15) and selectively including the evanescent states(we will 
use the right electrode as the example, the same is true for the left electrode): (1) We will 
start with all the running wave states, calculate the overlap matrix elements < ^^rj 
< (pn*\'Pm < 0n*l0m* ^^R projcctious < >Qii, then solvc the linear equations 



for A^{1), B^{1). (2) We will calculate the integral of wavefunction square of the right and 
left hand sides of Eq(15) within VLr, and calculate the percentage of the right hand side vs 
the left hand side results. We will call this decomposition percentage (which is alway less 
than 1). If this percentage is close to 1 within a criterion (e.g., 10~^), then stop. Otherwise 
go to next step. (3) We will include the next evanescent state which has its En{ke) closest 
to E. We will include (pn,e{r) = Un^k^expiikez) in summation of Eq(15), just treat it as one 
of the running wave states (but if ke is V or X', </>* e(r) is the same as (f)n,e{r), thus should 
not be included). Then repeat step (1),(2), find the new A^{1), B^{1), also the values for 
the evanescent states A^g(/), B^^{1). If one evanescent state has almost zero (e.g., less than 
10~^) contributions in all tp(^i){r), then discard this evanescent state. If the decomposition 
percentage is still not close enough to 1, repeat step (3). If the total number of evanescent 
state is too big (e.g, larger than 10), or the next closest En{ke) is too far away from E (e.g, 
farther than 2 eV), then stop. 

Let's assume that through the above procedure, we have included Nl, evanescent 
states (counting both possible 0„,e and 0* g) to the Ql and Qr sub-equations in Eq(15). If 
we have M > + Nl + Nf^ + (situation I), then we can request all the evanescent 
state coefficients to be zero after the linear combination in Eq(17) [e.g, Y^fii An^P{l)Ci = 0, 
Y.iLiB^'^^\l)Ci = ]. These are Nl + A^^ additional equations to Eq(18). Since we still 
have more number of Q than the total number of linear equations, we can still request 
the Y^f^ \Ci\^ /uji to be minimum while these equations are satisfied (again, this can be 
solved by the ZGELSS LAPACK routine [l^). Here we have placed a weight function 
uji, which depends on the decomposition percentage (after the inclusion of the evanescent 
states) of each ip{i) in Eq(15). If the decomposition percentage is close to 1 [a good fit 
in Eq(15)], then uji is close to 1. If the decomposition percentage is much less than 1 
[not a very good fit in Eq(15)], then uoi is very small, which means ■?/'(;) is discouraged 
from participating in the linear combination of Eq(16). More specifically, if we use and 
Pi to denote the decomposition percentage of the ip{i) at VLr and VLl in Eq(15), then we 
have used a formula ui = 0.001/(0.001 + bf - 1|) + 0.001/(0.001 + |pf - 1|). In another 
situation (situation II), we have Nr + Nl < M < Nr + Nl + Nf^ + N^ Then to solve 
Ci, we can minimize the evanescent state coefficients after the linear combination of Eq(16) 
[i.e, minimize En,H,L I Efii A^^^ {l)Ci\'' + \ EHi ^^i^HOCd'] while satisfying Eq(18) exactly. 
This again can be solved by standard numerical packages. In our calculation, we find all of 



our cases fall into situation 1. 

In the situation I, we have requested the evanescent state coefficients in the scattering 
state of Eq(16) to be zero. This might look strange at ffist. As we discussed above, the 
type I evanescent state might be physical in a scattering state. Then, how can we force 
it to be zero and still have a good scattering state? The answer lies in the fact that we 
did not separate the evanescence state (j)n,e{f') = Un,ke^xp{ikez) into the type I and type II 
states [ i.e., (pn,eG'^^ and </)„,ee~''^] . As a result, the coefficient we have for (t)n,e{,f) is really the 
sum of the coefficients for the type I and type II states. As a result, in a scattering state 
Ipse, although we cannot force the coefficients of type I evanescent states to be zero, we can 
always add a type II states to cancel their coefficients. So, when we require the coefficients of 
0n,e(^) to be zero, it doesn't mean the type I evanescent state coefficients are zero. However, 
since the possible type I evanescent state coefficient in a given scattering state ipse is fixed 
and is likely small sXVLl and VLr-, then our type II evanescent state coefficient should also be 
small. This guarantees that the erroneous situation of large type II evanescent states (as we 
discussed near the beginning of this section) will never happen, and our results are always 
stable and accurate. 



IV. THE RESULTS 



Following the above procedures, we have calculated the system in Fig.l with different 
biases. We have compared the current results with the results reported in Ref. 0. First, 
using the conventional ground state conjugate gradient program ,14], we have solved all the 
eigen states of H up to ~0.5 eV above the right electrode Fermi energy /i/j. In our system, 
this amounts to ~140 eigen states. Then we have scanned the scattering state energy E 
with an interval of ~ 0.04 eV. For each / and E, as shown in Fig.2, Eq(ll) of ip^^-^ can 
be solved by the conjugate gradient method within 20 line minimizations up to 10^^ a.u. 
accuracy. Next, decomposition of is carried out at VLr and VL^ as described by Eq(15) 
including the evanescent states. We find that, for most ■?/'(;), the running wave alone can 
get a decomposition percentage up to 99.999% or higher. However, for each energy E, it is 
very likely that there are one or two ?/'(/) with their running wave decomposition percentage 
only up to 50% or smaller. It is also likely that, even after including the approximated 
evanescent states, there are still one or two ipii) with their decomposition percentage only 



being around 90%. However, since these ipn) have very small ui in the minimization of 
Y,f^ ICip/u;;, their Ci are often exceedingly small (e.g, < 10^^°) in the linear combination of 
Eq(15). Following the procedure described above, a scattering state ipse of Eqs(16),(17),(18) 
is solved for each right electrode running wave 0^*. The typical ipse wavefunctions look the 
same as illustrated in Fig. 3 of Ref. jsj. The transmission coefficient Tm{E) is calculated for 
the scattering state according to Eq(4). We find that the current conservation equation (5) 
is mostly satisfied beyond 99.9%, and in many cases beyond 99.999%, an indication of the 
numerical accuracy of this approach. However, there are occasional and distinctive cases 
where Eq(5) is not satisfied at all (e.g., the sum of transmission and refiection is 10^, instead 
of 1). In these cases, the evanescent states in the constructed scattering state of Eq(16) are 
not eliminated, but dominating, perhaps due to our approximations in the treatment of the 
evanescent states. Fortunately, these cases are very rare and can be easily detected, thus to 
be discarded. 

Figure 4 shows the calculated transmission coefficients Tn{E) for different band n, plotted 
as functions of the kz points. The system has a bias of 4V. Each cross symbol corresponds 
to a calculated scattering state ipse- We have also plotted the calculated Tn{E) using our 
previous method 8] as rectangular symbols. As we can see, the current method and the 
previous method yield the same T„ amplitudes. This is a cross check for the robustness of 
these two methods. However, since the E in the previous method can only be the eigen 
energies Ei of the original H, it only yields a finite number of the scattering states. In 
contrast, the current method can have as many scattering states as we want. Actually, as 
shown in Fig.4, there are cases (e.g, for n = 5 and for the dip near kz = O.dir/a of n = 1) 
where the previous method might not have enough calculated points to reveal the sillenXX 
nature of the Tn{kz) curve. In terms of the computational time, we find the current method 
is faster than the previous method, despite the fact that we now have much more data points. 
We find that for the number of E points we used, the total time spent to solve Eq(ll) for 
all the energies and I is about 2 times the time spent to solve all the ground states ipi of the 
original Hamiltonian H. This makes our transport computational time in the same order as 
a typical ground state calculation. 

The points in Fig.4 are fitted by smooth curves Tn{kz) as described in Ref. 0, and the 
resulting curves are used to plot the total transmission coefficients T{E) = J2nTn{kz{E)), 
which is shown in Fig. 5. Again, we have compared our current results with our previous 



results j8[ for the biases of 1 V and 4 V cases. We see that, overall, they are almost the 
same. But there are some differences in the detail. In the bias 1 V case, near E — Ep = —1.5 
eV, the current method produce a well shape dip. This is due to a gap at the X' point of 
the left electrode near E — fi^ = —0.5 eV. The previous method missed this dip because it 
doesn't have a data point with its energy falls into this left electrode energy gap. In the 
bias 4 V case, near E — Ep = eV, the current T{E) is lower than the previous results. 
This is because in the n = 5 band of Fig. 4, the previous method has only two points at the 
kz < O.Gvr/a region. This leads to a fitted Tn{kz) which is too high compared to the correct 
result, and consequently an over estimated T{E) near E — Ep = 0. 

Despite the above differences between the current T{E) and the previous results, their 
calculated total currents are very similar. For example, in the cases of 1 V and 4 V biases, the 
current method produces currents 0.0390 and 0.376 e^V/h respectively, while the previous 
method produces 0.0417 and 0.398 e^V/h respectively. The differences are only about 5%. 
The I-V curve produced by the current_jiiethod is very close to the result of the previous 
method, which is shown in Fig. 6 of Ref 



V. CONCLUSION 



In conclusion, we have presented a new approach to calculate the quantum transports. 
The current method is based on a previous method Q which uses the periodic boundary 
conditions, thus makes it possible to use the popular pseudopotentials and planewave basis 
set. Compared to the previous method however, the current method uses a different way 
to solve the periodic wavefunction ■?/'(/) of Eq(6). As a result, the scattering states can be 
calculated at any given energy E. This provides a more robust way to calculate the scattering 
state wavefunctions and their transmission coefficients. Under the current method, the total 
computational time for a transport problem is in the same order as the computational time 
of its corresponding ground state problem. Enough details of the procedure is presented 
which makes the implementation of this method possible. 
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FIG. 1: A schematic view of the calculated system. 



FIG. 2: The conjugate gradient (CG) convergence of Eq(ll). The convergence error is defined as 
II ||. 

FIG. 3: The band alignment between the left electrode band structure and the right electrode 
band structure. The voltage bias is 4 V. The numbers in the right electrode band structure are the 
band index. The small black dots on the line E are the and points satisfying E^{k^) = E 
and E^{k^) = E respectively. 



FIG. 4: The calculated transmission coefficients Tn{k^). The crosses are the results from the 
current method, the rectangulars are the results from the previous method Ref. 0, and the lines 
are the fitted smooth curves for the current results. 



FIG. 5: The calculated total transmission coefficients T(E) for different biases. The zero is the 
right electrode Fermi energy. For a given bias V, there are net right to left current flow only within 
the [— V^,0] energy window. 
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